Electron correlation in C4/V+2 carbon rings: aromatic vs. dimerized structures 
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The electronic structure of Cijv+2 carbon rings exhibits competing many-body effects of Hiickel 
aromaticity, second-order Jahn- Teller and Peierls instability at large sizes. This leads to possible 
ground state structures with aromatic, bond angle or bond length alternated geometry. Highly 
accurate quantum Monte Carlo results indicate the existence of a crossover between Cio and C14 
from bond angle to bond length alternation. The aromatic isomer is always a transition state. The 
driving mechanism is the second-order Jahn- Teller effect which keeps the gap open at all sizes. 



The discovery of carbon fullerenes and nanotubes has 
opened a new materials research area with a vast poten- 
tial for practical applications. Unfortunately, our under- 
standing of the rich variety of structural and electronic 
properties of carbon nanostructures is far from satisfac- 
tory. For example, experiments B indicate that quasi- 
one-dimensional structures such as chains and rings are 
among the primary precursors in the formation process 
of fullerenes and nanotubes. However, our insights into 
their properties and behavior are incomplete due to the 
complicated many-body effects involved, fn particular, 
recent studies [^]|| have demonstrated a profound im- 
pact of the electron correlation on stability and other 
properties of such all-carbon structures. An important 
example of such nanostructures is the system of planar 
monocyclic carbon rings C„ with n = AN + 2, where N 
is a natural number. These closed-shell molecules mani- 
fest an intriguing competition between conjugated aro- 
maticity, second-order Jahn- Teller and, at large sizes, 
Peierls instability effects. Consequently, this leads to dif- 
ferent stabilization mechanisms that tend to favor one of 
the following structures: a cumulcnic ring (A), with full 
T) n h symmetry, with all bond angles and bond lengths 
equal; or either of two distorted ring structures, of lower 
D.jh symmetry, with alternating bond angles (B) or bond 
lengths (C). Further structural details are given in Fig. [j]. 
Accurate studies for the smallest sizes (Cg and Cm) find 
isomer B to be the most stable. However, for larger sizes 
the results from commonly used methods are contradic- 
tory and available experiments Q are unable to clearly 
specify the lowest energy structures. 

In order to identify the most stable isomers and to elu- 
cidate the impact of many-body effects, we carried out 
an extensive study of electronic structure and geometries 
of C4AT+2 rings of intermediate sizes up to 22 atoms (with 
some methods up to 90 atoms). We employed a number 
of electronic structure methods including the highly ac- 
curate quantum Monte Carlo (QMC) method which has 
been proven very effective in investigations of C20 and 
larger carbon clusters ||, as confirmed also by an inde- 
pendent study by Murphy and Friesner |J . Our QMC re- 
sults reveal that the ground state structures have 
alternated geometries at all sizes while cumulcnic isomer 



A is a structural transition state. The results also pro- 
vide valuable insights into the shortcomings of the den- 
sity functional approaches such as inaccurate balance be- 
tween exchange and correlation in commonly used func- 
tional. In addition, the letter presents a first evaluation 
of interatomic forces in large systems within the QMC 
framework. 

According to the Hiickel rule, the 4A+2 stoichiometry 
implies the existence of a double conjugated 7r-electron 
system (in- and out-of-plane). Combined with the ring 
planarity, this suggests a strong aromatic stabilization in 
favor of isomer A. Although the highest occupied and 
the lowest unoccupied molecular orbitals (HOMO and 
LUMO) are separated by a gap of several eV, a double 
degeneracy in the HOMO and LUMO states opens the 
possibility for a second-order Jahn- Teller distortion [n|, 
resulting in either cumulenic B or acetylenic C structure. 
Such distortion lowers the symmetry and splits the de- 
generacy by a fraction of an eV, with an overall energy 
gain. Moreover, as N — > 00, the picture is complicated 
further by the fact that the system becomes a semimetal- 
lic polymer with two half-filled it bands. As first pointed 
out by Peierls , such a one-dimensional system is intrin- 
sically unstable and undergoes a spontaneous distortion 
which lowers the symmetry. The symmetry breaking en- 
ables the formation of a gap, in analogy to the elusive 
case of irems-polyacetylene § . 




FIG. 1. The most stable isomers of C4JV+2 rings (shown 
for Cio). The parameters needed to specify the geometries 
are: average bond length r = (n + ?"2)/2, and bond angle 
a — (qi + Q2)/2; relative bond length A r = (ri — r2)/f, and 
bond angle alternation A a = (ai — 012) /a. 
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It is very instructive to see how the commonly used 
computational methods deal with such many-body ef- 
fects. Density functional theory (DFT) methods tend to 
favor a "homogenized" electronic structure with delocal- 
ized electrons. In fact, for sizes larger than Cio, there is 
no indication of any stable alternation up to the largest 
sizes we have investigated (C90). Calculations performed 
within the local density approximation (LDA) and gen- 
eralized gradient approximations (GGA, with BPW91 
functional) consistently converge to the aromatic struc- 
ture A, in agreement with other studies S. Only by 
extrapolation to the infinite-chain limit, Bylaska, Weare 
et al. jlo| claim to observe a very small, yet stable, bond 
alternation within LDA. A very different picture arises 
from the Hartree-Fock (HF) method, which shows a pro- 
nounced dimerization for Cio and larger. This agrees 
with the HF tendency to render structures less homoge- 
neous in order to increase the impact of exchange effects. 
We also verified that using GGA functionals with an ad- 
mixture of the exact HF exchange (B3PW91) recovers 
qualitatively the HF results for large sizes (>C46), as al- 
ready observed by others . 

Obviously, this problem calls for much more accurate 
treatments. High-level post-HF methods, such as multi- 
configuration self-consistent field (MCSCF) and coupled 
cluster (CC), indeed provide answers for the smallest ring 
sizes (C 6 pi and Cm @@)- In particular, Martin 
and Taylor [jl4[ have carried out a detailed CC study 
demonstrating that both C6 and Cio have angle alter- 
nated ground state structures, although for Cio the en- 
ergy of the aromatic isomer A is found to be extremely 
close (1 kcal/mol). In addition, we have performed lim- 
ited CCSD calculations of C14 and have found the dimer- 
ized isomer to be stable by ~6 kcal/mol. Unfortunately, 
these methods are impractical for larger cases or more 
extensive basis sets [ pT| . 

The quantum Monte Carlo (QMC) method was used 
to overcome these limitations. This method possesses 
the unique ability to describe the electron correlation ex- 
plicitly and its favorable scaling in the number of parti- 
cles enables us to apply it to larger systems fig] . In the 
variational Monte Carlo (VMC) method we construct an 
optimized correlated many-body trial wavefunction 
given by the product of a linear combination of Slater 
determinants and a correlation factor 



n I,i<j 



u(rn,rji,rij) (1) 



where tp are one-particle orbitals, i,j denote the elec- 
trons, / the ions and Tu,rji,rij are the corresponding 
distances. The correlation part, it, includes two-body 
(electron-electron) and three-body (electron-electron- 
ion) correlation terms and its expansion coefficients are 
optimized variationally. Most of the variational bias 
is subsequently removed by the diffusion Monte Carlo 
(DMC) method, based on the action of the projection 
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FIG. 2. HF and MCSCF energy as a function of the degree 
of dimerization A r - Least-squares estimates of the positions 
of the minima are indicated by arrows. 

operator exp(— tH); in the limit of r-^oo, this projec- 
tor recovers the lowest eigenstate from an arbitrary trial 
function of the same symmetry and nonzero overlap. The 
fcrmion antisymmetry (or sign) problem is circumvented 
by the fixed-node approximation. More details about the 
method are given elsewhere |T^] . 

DFT, HF and MCSCF calculations have been carried 



out using standard quantum chemistry packages 16 1. All 
calculations employed an accurate basis set, consisting of 
10sllp2d Gaussians contracted to 3s3p2d, and smooth ef- 
fective core potentials {l?]] to replace the chemically inert 
core electrons. 

The geometries of smaller rings with 6 and 10 atoms 
have already been established from previous calcula- 
tions |L2-14|. We have verified that the most reliable 
published structural parameters agree very well (within 
~ 0.002 A and 1°) with our own GGA values. However, 
since the dimerized isomer C is unstable within DFT, 
we followed a different strategy. We began from HF ge- 
ometries, which show that the degree of bond length al- 
ternation saturates at A r ~ 14% (Fig. ||). In order to 
correct for the HF bias favoring acetylenic structures, we 
performed limited MCSCF calculations (see below) for 
Cio, C14, and Cis- The electron correlation has a pro- 
found effect on the geometry, to the extent of causing 
the dimerized isomer to be unstable for Cio, while for 
C14 it decreases the dimerization to A r ~ 10%. Clearly 
the limited MCSCF for C14 and Cis provides rather poor 
geometry improvement although one expects a larger cor- 
rection as more correlation energy is recovered. In order 
to verify this and to estimate the correct degree of dimer- 
ization for C14, we carried out the evaluation of the Born- 
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Oppenheimer forces by a finite-difference scheme using 
correlated sampling, in the VMC method Jl^Jl^]. The 
computation of interatomic forces is a new development 
in QMC methodology and, to our knowledge, this is the 
first application in this range of system sizes. We probed 
the tangential C-C stretching/shortening which leads to 
a change in the degree of dimerization, A r - For A r = 7%, 
our force estimate is F — —dE/dO = 0.010(2) a.u. (and 
a second derivative of H = 0.30(1) a.u.), suggesting 
proximity to the minimum. Moreover, at A r — 10.5% 
we find a force of opposite sign: F = —0.013(3) a.u. 
(H = 0.33(1) a.u.). For Cis, we have instead performed 
two QMC single point calculations at A r — 7%, 14% and 
found the first energy to be lower by AE ~ — 12 kcal/mol. 
Finally, we assumed A r — 7% and r = 1.286 A as our best 
estimate for calculations of the acetylenic isomer with 
n > 10. 

The crucial ingredient for very accurate QMC calcu- 
lations is a trial function with a small fixed-node error. 
The quality of the one-particle orbitals is of prime im- 
portance for decreasing such error. While HF or DFT 
orbitals are commonly used for construction of Slater de- 
terminants, our recent projects p(^] have demonstrated 
that natural orbitals from limited correlated calculations 
(e.g., MCSCF) lead to more consistent results. Inclu- 
sion of the electron correlation into the method used 
to generate the orbitals is indeed very relevant for ob- 
taining improved fcrmion nodes, especially for such sys- 
tems which exhibit strong non-dynamical correlation ef- 
fects |y],[l4| • Extensive tests confirmed that orbitals from 
MCSCF (with an active space consisting of 4 occupied 
and 4 virtual orbitals) yield the lowest energies and so 
they were used in all of our calculations. In addition, 
the inclusion of the most important excited configura- 
tions into (about 20-30 determinants) provided fur- 
ther significant improvement of the total energies. In 
particular, the weights of single excitations were surpris- 
ingly large for the alternated geometries and comparable 
to the largest weights of configurations with double exci- 
tations. A more detailed analysis on the multi- reference 
nature of the wavefunction in these systems will be given 
elsewhere. 

Equipped with such high quality trial functions we 
have carried out QMC calculations from Gq to Cig. A 
plot of the energy differences, with comparison to other 
methods, is shown in Fig. ||. For the very compact C§ 
ring, where the overlap between 7r in-plane orbitals is 
large, as observed by Raghavachari et al. pl| , the angle 
alternated isomer B is the most stable. The aromatic 
structure A is instead a transition state leading to angle 
alternation (Bi„ mode), while the dimerized isomer C is 
unstable in all methods. 

Cio is the case which was studied extensively in the 
past. Our DMC results agree with calculations of Martin 
and Taylor We conclude that the angle alternated 
isomer is still the lowest energy structure, albeit 
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FIG. 3. Relative energies of the angle alternated B (filled 
symbols) and bond alternated C (empty symbols) structures 
with reference to the aromatic isomer A (see Fig. Q). 

extremely close to the cumulenic A geometry, with an 
energy difference of ~ 1 kcal/mol. Indeed, the stabiliza- 
tion by aromaticity is almost as strong as the effect of 
second-order Jahn- Teller distortion which is responsible 
for the alternation pattern. The aromatic isomer remains 
a transition state, as it is for C6, although in this case 
the energy surface is extremely flat. The Cio acetylenic 
isomer appears unstable in DMC, which implies that our 
older calculations j2Q] and also a more recent all-electron 
study p2[ |, based on a single determinant vpy with HF 
orbitals, were not accurate enough. 

Our perhaps most interesting results come from the 
C14 and Cis isomers. The angle alternated structures be- 
come unstable since the in-plane orbital overlap is smaller 
due to the increased ring radius. The trends from HF, 
MCSCF and QMC (Fig. |) are clearly in favor of dimer- 
ized geometries, although there is an indication that the 
HF bias for the Jahn- Teller distortion is much reduced 
as we recover more of the correlation energy. Neverthe- 
less, since the fixed-node DMC method typically obtains 
~ 95% of the correlation energy, we argue that margins 
for possible errors are very small. This is in sharp con- 
trast with the density functional results, which indicate 
only the aromatic isomer A to be stable. It seems that 
DFT methods "overshoot" the effect of correlation at the 
expense of exchange, resulting in a qualitative disagree- 
ment with QMC results. 

The data from HF and QMC calculations enable us to 
analyze the HF and correlation energy components sep- 
arately as a function of ring size. For n > 10, the HF 
energy difference between aromatic (A) and dimerized 
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(C, A r = 7%) isomers can be approximated by 

e hf - e hf - 5-14n - 45.4 [kcal/mol]. (2) 

For the correlation energy such an extrapolation is less 
certain, as we lack data for very large sizes. The fol- 
lowing formula reproduced our C10-C22 data within the 
error bars obtained (the value for C22 was based on an 
adjusted single determinant DMC calculation) 

Ecorr - E^rr - -3.57n + 27.4 [kcal/mol] (3) 

while for larger sizes it is probably somewhat less ac- 
curate. Nevertheless, the significant difference between 
the coefficients enables us to recognize the dominance 
of the Coulomb and exchange contributions for large n. 
This means that the competition between the stabiliza- 
tion mechanisms, which determine the geometry, gap, 
and key electronic structure features, is effectively de- 
cided at intermediate ring sizes. On the basis of these 
results we suggest that the ground state isomers are al- 
ternated at all sizes. 

In conclusion, we present a study of C47V+2 carbon 
rings, for which highly accurate fixed-node diffusion 
Monte Carlo results show that distorted isomers are al- 
ways preferred with a stability crossover from bond angle 
to bond length alternation between C10 and C14. The 
fully-symmetric aromatic isomer is instead a structural 
transition state connecting equivalent alternated geome- 
tries shifted by one site along the ring, although for C10 
the angle alternated isomer is below the aromatic one by 
only ~ 1 kcal/mol. The intermediate size rings (above 
C10) show a clear trend of dimerization, indicating that 
the driving stabilization mechanism is the second-order 
Jahn- Teller effect. The corresponding HOMO-LUMO 
gap persists in large rings and overshadows thus the on- 
set of the Peierls regime. The calculations also provide 
an interesting insight into the deficiencies of density func- 
tional approaches such as the imbalance in approximat- 
ing the exchange and correlation energy contributions. 
The competition of complicated many-body effects, from 
aromatic conjugation to symmetry-breaking, provides an 
extremely sensitive probe of the treatment of electron 
correlation in different theoretical approaches. Our re- 
sults demonstrate the potential of quantum Monte Carlo 
not only for evaluation of accurate energy differences, but 
also for other properties, such as equilibrium structures, 
which can be obtained by new developments in the QMC 
methodology. 
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